BioVCN R Lesson 07 – Statistics

Today we are going to talk about running some basic statistics in R.

library(tidyverse)
library(plotly)
library(cowplot)

********************************************************
Note: As of version 1.0.0, cowplot does not change the
  default ggplot2 theme anymore. To recover the previous
  behavior, execute:
  theme_set(theme_cowplot())
********************************************************

Usually, when I want to run stats on a thing, I start with some basic descriptive stats. Fortunately we did this last week with the starwars dataset, so lets jump back in.

data(starwars)

Lots of statistics are applied to continuous variables. Lets take a look at the relationship between height and mass. I’m going to make an interactive plot with ggplotly so we can mouse over points and see them better.

Here I included a name and spec variable, but didn’t map them to anything. Thats so you can see them when you mouse over them.

hmPlot <- starwars %>% ggplot(aes(x = height, y = mass, name = name, spec = species)) + geom_point() + theme_cowplot()
ggplotly(hmPlot)

Woah dude. One character is anomalously heavy. Mouse over the outlier so you can see who it is.

That was my first guess too

Jabba is probably going to screw up all of our stats, lets focus on the relationship of all of the other charactersby removing him.

starwars01 <- starwars %>% filter(str_detect(name, "Jabba")) # Complecated syntac because I can't spell the rest of his name

hmPlot01 <- starwars01 %>% ggplot(aes(x = height, y = mass, name = name, spec = species)) + geom_point() + theme_cowplot()
ggplotly(hmPlot01)

So here’s our relationship. It looks sort of, but not reeally linear, which makes sense, there are lots of species in the galaxy.

Sometimes when we do stats, the stats like to imagine normally distributed data.

starwars01 %>% pivot_longer(cols = mass:height, names_to = "measurement", values_to = "value") %>% 
  ggplot(aes(x = value)) + facet_wrap(measurement~.) + geom_histogram()

Hmm. Sort of.

See if you can figure out what I did above.

Lets see if height and width are correlated.

cor.test(starwars01$height, starwars01$mass, method = "pearson")

    Pearson's product-moment correlation

data:  starwars01$height and starwars01$mass
t = 8.7853, df = 56, p-value = 4.018e-12
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.6260700 0.8520232
sample estimates:
      cor 
0.7612612 

cor.test(starwars01$height, starwars01$mass, method = "spearman")

Cannot compute exact p-value with ties

    Spearman's rank correlation rho

data:  starwars01$height and starwars01$mass
S = 6946.8, p-value = 2.597e-13
alternative hypothesis: true rho is not equal to 0
sample estimates:
     rho 
0.786313 

Here are a parametric and non-parametric correlation test. The values I usually think about are cor (pearson), which is the R value, or rho(spearman) which is its “rho” which is a lot like an R value.

Numbers closer to 1 are more positively correlated, closer to -1 are more negatively correlated, closer to 0 are not correlated.

There is a p-value for both of these, which is essentially the probability that if all of the assumptions of the model hold, the true R or rho value is zero.

Regressions

Regressions are also pretty popular with the kids these days. Lets do one! Lets see how well height explains mass

mod <- lm(mass ~ height, data = starwars01)
summary(mod)

Call:
lm(formula = mass ~ height, data = starwars01)

Residuals:
    Min      1Q  Median      3Q     Max 
-39.382  -8.212   0.211   3.846  57.327 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -32.54076   12.56053  -2.591   0.0122 *  
height        0.62136    0.07073   8.785 4.02e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 19.14 on 56 degrees of freedom
  (28 observations deleted due to missingness)
Multiple R-squared:  0.5795,    Adjusted R-squared:  0.572 
F-statistic: 77.18 on 1 and 56 DF,  p-value: 4.018e-12

Lets look at what this tells us. Residuals are essentially the distance of each point from the models prediction. The one that is farthest below is ~ -39 too low, the highest above is ~ 57 and the median and first and third quartiles are also shown.

Coefficients are often what you care about. The estemate tells us about the slope and intercept of the linear regression. These values basically say our model looks like this

mass = -32.5 + 0.62 * height

We don’t know the true values of these paramters, so the standard error gives us an idea of the ranges of the two. T-value is the test score to see how good those are, band the p-value tells you if the answer is statistically significant.

4.02 x 10^-12 < 0.002 so it the height value gets a bunch of stars after it. I usually don’t make too much of the p value of the intercept.

Residual standard error is the standard deviation of the residuals, according to the internet. 56 degrees of freedom is calculated as our sample size, minus the complexity of the model. We have 58 characters, but we loose one DF for the intercept and another for the height varaible that we are using.

There is also an R^2, and adjusted R^2 which gets smaller as the model gets more complicated, and a p-value for the whole model.

Residuals

Linear models like to assume that the residuals are normally distributed. We want to make sure that no person is really affecting our model too much. We can plot the model to do this.

Processing models

the broom package gives a nice matrix version of model results

library(broom)
tidy(mod)